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Media composed of colliding hard disks (2D) or hard spheres (3D) serve as good 
approximations for the collective hydrodynamic description of gases, liquids and 



bX)' granular media. In the present study, the compressible hydrodynamics and shock 

<^ . dynamics are studied for a two-dimensional hard-disk medium at both the contin- 

uum and discrete particle level descriptions. For the continuum description, closed 
form analytical expressions for the inviscid hydrodynamic description, shock Hugo- 
^ . niot, isentropic exponent and shock jump conditions were obtained using the Helfand 

equation of state. The closed-form analytical solutions permitted us to gain physi- 
^ ■ cal insight on the role of the material's density on its compressibility, i.e. how the 

medium compresses under mechanical loadings and sustains wave motion. Further- 

>v 

more, the predictions were found in excellent agreement with calculations using the 
!— Event Driven Molecular Dynamic method involving 30,000 particles over the en- 
tire range of compressibility spanning the dilute ideal gas and liquid phases. In all 
cases, it was found that the energy imparted by the piston motion to the thermal- 
ized medium behind the propagating shock was quasi-independent of the medium's 
packing fraction, with a correction vanishing with increasing shock Mach numbers. 
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I. INTRODUCTION 



Dilute media composed of colliding hard disks (2D) and hard spheres (3D) serve as 
very good molecular models for gases, liquids and granular media. Indeed, the kinetic 
theory of dilute hard sphere gases is well established^, whereby a hydrodynamic or coarse- 
grained description, such as the Navier-Stokes equations, can be obtained from the collective 
dynamics of many particles. When the density of the hard particle media is increased, such 
that the distance travelled by a particle before experiencing another collision (i.e., the mean 
free path) becomes comparable with the dimension of the particle itself, departures from the 
ideal gas behavior become important. For such regimes, hard particle media more closely 
resemble liquids and even solids. Significant work has thus been devoted to determine the 
equation of state of dense hard particle medial, phase transition phenomena^- and shock 
wave structure^ - —. With the advent of more powerful computers, hydrodynamic problems, 
such as the shock wave implosion,— have been studied using hard sphere models. Likewise, 
the hard particle description also permits to study complex problems in chemical kinetics 
and reactive dynamics^— including detonation wave propagation^^. In all cases, the hard 
particle paradigm significantly simplifies the collision dynamics, such that kinetic-theory and 
statistical mechanics tools are more easily implemented to pass from the molecular (particle) 
dynamic description to the continuum description at larger scales. 

One of the recent successes of the hard particle dynamics paradigm is to demonstrate 
how large hydrodynamic fluctuations can arise in granular media if one takes the collisions 
to be endothermic (dissipative)^. Indeed, large scale clustering, spontaneous vortical flow 
formation, super- diffusion and departures from Maxwell-Boltzmann distributions arise when 
the collisions among particles are inelastic. These findings have triggered a renewed interest 
in the kinetic theory of hard particle media; a significant literature is now available on the 
subject, recently reviewed in Refs. 18 and 19 . The framework, however, is usually based on 
the dilute limit, although granular flows approach the structure of a liquid, or even solid, 
for which the compressibility is likely to severely affect the predictions. The present study 
addresses the dense regime of hard particle hydrodynamics. 

Granular flows of macroscopic particles are typically highly compressible, in that the 
granular sound speed, which is approximately the mean speed of the particles, is usually 
comparable or smaller than the velocity of the solid surfaces driving the system. For exam- 
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pie, consider the walls of a vibrated container— ~—, whose walls impart their energy to the 
colliding particles, or the propeller blades of a mechanical mill agitating grinding balls^. 
Likewise, the chute and avalanches of granular particles over obstacles 2 ^ 1 ^ may also generate 
strong shock waves, while the granular temperature and sound speed of the particles is very 
low owing to the low thermal velocity of the particles. Such flows are thus characterized 
by the propagation of strong shock waves, propagating at large Mach numbers, often in the 
hypersonic regime. The flow regime of granular media is thus highly compressible. 

The present study focuses on the compressible dynamics of dense hard particle media, i.e. 
media in which the ideal gas assumption fails. We focus on the description of inviscid flows in 
such media, with a particular emphasis on shock wave solutions. Shock wave propagation, as 
generated by the sudden motion of a piston into a thermalized medium, is directly analogous 
to the sudden acceleration of a wall in a granular system^ or to the classical problem of a 
shock generated by a moving piston in gas dynamics^. We focus on large departures from the 
dilute limit and wish to describe the energy addition to the system by the piston's motion. 
Because we wish to compare the analytical results with molecular dynamic simulations 
and wish to easily visualize the system's dynamics, the present work focuses on the hard- 
disk system in two- dimensions. This also permits to conduct calculations on much larger 
scales. The restriction of our study to two space dimensions does not however affect the 
generality of the method and results, as they can be simply extended to three dimensions 
with the appropriate equation of state. The present study also assumes the collisions to be 
elastic. The shock propagation problem through a dissipative medium, for which continuum 
descriptions have been proposed in Ref.— , is left for future study. Instead, we focus first 
on developing and validating the use of a sufficiently simple and accurate equation of state 
given by Helfand et al.— , for which analytical solutions can be obtained for the dense regime 
offering a clear insight into the role of material's density on compressibility. This will permit 
future extension to the dissipative regime. 

The paper has the following structure. In the first part, we use the simple equation of 
state for a hard disk fluid proposed by Helfand and Frisch 28 to develop analytical descriptions 
for the medium's inviscid hydrodynamics and shock wave solutions. In the second part of 
the paper, we report the results of molecular dynamic simulations and compare with the 
analytical predictions, hence further validating the model derived and further clarifying the 
role of the medium's density on its compressibility and shock dynamics. 
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II. CONTINUUM DESCRIPTION 



A. Equation of state for a hard disk gas 

In a hard particle system, the internal energy e of the system consists only of the trans- 
lational kinetic energy modes, with |— T per translational degree of freedom, where k is the 
Boltzmann constant, m the particle mass, and T the temperature. For a hard disk system, 
there are only two translational degrees of freedom, yielding the caloric equation of state 
given by 

An equation of state relation for the thermodynamic state variables (p, v, T) also needs to 
be prescribed; generally, it can be written as^: 

pv = Z—T (2) 
m 

where Z is the compressibility factor, which is unity for an ideal (dilute) gas. Combining 
(PQ) and (j2J), the internal energy can be written uniquely in terms of p and v. 

e = Y (3) 

The compressibility factor Z is usually expressed in terms of the packing factor 77, denoting 
the fraction of the volume occupied by the particles, i.e. 

Va/m 

V = (4) 

V 

where V a is the volume (i.e. the surface area of the disk in 2D) of a particle. The caloric 
equation of state is thus uniquely expressed in terms of the pressure and specific volume in 
the form 

pv 

eM = m^j (5) 

For a hard-disk medium, an exact equation of state is not known in closed form, hence 
requiring an infinite virial expansion. A thorough discussion of the various equations of 
state proposed for a hard disk system and their merits are described in Ref.-. An equation 
of state that is both accurate, physically meaningful and analytically simple is that proposed 
by Helfand et al.— . It takes the form 




FIG. 1. The relationship of the compressibility factor, with the packing factor 77 

To highlight its accuracy, we have compared it with the more accurate equation of state 
proposed by Maeso et alM, which takes the form 

10 3 - 947.989?? + 128.018r7 2 - 113.98777 3 



Z( V ) 



;i- 0.947989??) (1-77) 2 10 3 
-52.9722r? 4 - 1.58059677 5 



(l- 0.947989r/)(l-?7) 2 10 3 1 ' 

This equation of state is obtained from the first seven virial coefficients and has an accuracy 
characterized by an average absolute deviation of 0.05%. The two equations of state are 
compared in Fig. [TJ As can be verified from Fig. [TJ Helfand's equation of state offers a 
very good compromise between simplicity and accuracy over the entire range of compaction 
factors 77. It is generally accurate to within ~ 3% across the entire range of packing fractions 
spanning the dilute gas, dense gas, liquid and solid states close to the crystal packing fraction 
corresponding to a regular triangular lattice, given by 

nV3 



Vc 



6 



0.907 



(8) 



We will henceforth adopt this equation of state for analytical simplicity and insight into 
the effect of initial density (packing factor) on the medium's compressibility. 



B. Isentropic exponent and sound speed 



The change in the compressibility of the medium is best described by the isentropic 
exponent 7, which describes the relation between changes in density and changes in pressure 
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for isentropic processes. It is defined as 

where p — 1/v is the density. The sound speed in the medium is directly related to 7, i.e.: 

2 _ ( d P\ . P 



^WrS (10) 

Note that 7 is not the ratio of specific heats, which in this case is simply 2 for any 
compaction ratio, as can be deduced from (JTJ and the definition of enthalpy. To determine 
the dependence of c or 7 on the state variables p and v (or p), we have to make use of the 
Gibbs thermodynamic equations relating p, v and the entropy s. The first and second law 
of thermodynamics for a control system yield Gibbs equation written in terms of the density 
as 

deOp, p) = Tds + —dp (11) 
p 2 

alternatively, e(p, p) can also be written as a perfect differential as 

deiP ' p)= {^) p dp+ {%) p dp (12) 



Comparing (ITT]) and (fT2"|) . we can re- write 

P _ ( de\ 



dp J \dp 

Now writing p(p, s) as a perfect differential, we get 



P 2 \dpj T 
dP= /ge\ P d P+7£T d8 ^ 



dp J s \ds 
Comparing (fT3"|) with (jHJ), we get immediately that 



dp(p, *)=[?£) dp+(^ \ ds (11) 



JP / de \ 



^'.Sj.-'TSr <15) 



P 



dp J 

This expression is general for any medium, as we have not yet used the equation of state of the 
medium. Using expressions (BJ, (jSJ) and (jBJ) to evaluate the derivatives, we get immediately 
the expression for the sound speed in terms of the local state of the medium and the local 
packing factor rj 

c 2 = pv (1 + (1 - p)- 2 + 277(1 - p)- 1 ) (16) 
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FIG. 2. The variation of the isentropic exponent with the local packing factor r\ 
and from ffTU]) we get the isentropic exponent. 

7 = 1 + (1 - r])- 2 + 2r7(l - J])' 1 (17) 

Figure [2] shows the variation of the isentropic exponent with the local packing fraction. In the 
dilute limit of an ideal gas, i.e. rj — > 0, we recover the isentropic exponent of a 2-dimensional 
hard disk gas of 7 = 2. However, with increasing packing fraction, the isentropic exponent 
grows commensurably, reflecting the incompressibility of the medium. A small change in 
density requires a very large change in pressure, a characteristic of nearly incompressible 
media like liquids and solids. 

We stress out that having written the isentropic exponent in closed form permits us to 
formulate the hydrodynamic equations for the medium very simply. This is performed in 
the next section. 



C. Gasdynamic description of a hard disk medium 

The compressible motion of a medium, after neglecting slow processes involving molec- 
ular transport, are governed by the Euler equations. Since the resulting motion is particle 
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isentropic, we can write them as: 



^ = -PV • u (18a) 

= ~ Vp (18b) 

Ds . . 

— = 18c 
Dt V ; 

where u is the macroscopic velocity vector and D/Dt = d t + u • V is the usual material 
derivative, following an element of fluid. Using (|14p . applying it on the particle path (i.e., 
following a thermodynamic system) and making use of f lTUj) and (Q, one can re-express f)18c[) 
as 

^M = 7W £^) (19) 



Using the relation for 7 given by (ITT)) , this can be integrated in closed form in order to 
determine the isentropic relation between pressure and density along a particle path. Taking 
the reference state as p\ and p%, the isentrope becomes 



P _ ( P \ ( pi 11 \ e (i-^i)d-,i) (20) 



Pi VPi/ \ 1 - Vi 

where we recover the result obtained by Gaspard and Lutskc 11 using a more complex equa- 
tion of state. For ideal hard-disk gases, in the limit of vanishing packing factor, we recover 
the isentropic relation p oc p 2 . 



D. Shock Waves in a Hard-Disk Medium 

When shock waves are generated in a medium, weak solutions to the inviscid Euler 
equations are required. Consider a piston suddenly accelerated from zero velocity to a 
constant velocity u = u p into a medium initially at rest; see Fig. [31 A shock wave is formed, 
which propagates with velocity D. If diffusive fluxes are negligible in the pre- and post-shock 
states, and letting the subscript 2 denote the uniform state of the medium behind the shock 
wave, the conservation of mass, linear momentum and energy, in the frame of reference of 
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piston shock 



D-u 



FIG. 3. Propagation of a piston driven shock wave into a quiescent medium and transformation 
to the shock fixed frame of reference (below) 



the shock (see Fig. [3]) yield 26 - 



D _ (D- u p 



D 2 

Pi H = P2 + 



V2 

(D ~ u p ) 2 
v 2 



ei + PiV! + ^D 2 = e 2 + P2V2 + -{D- u p ) 2 



(21) 
(22) 
(23) 



where v is the specific volume, p is the hydrodynamic pressure and e is the specific internal 
energy. Combining the equation of mass f )2ip . momentum (122 j) and energy (12 3 p to eliminate 
the two speeds, the accessible end states 2 across the shock wave can be represented by the 
shock Hugoniot relation in the (p, v) plane 2 ^: 



e2i.P2.v2) - ei(p 1( ui) = ~(vi - v 2 ){pi +p 2 ) 



(24) 



with either D or u p as parameters from (|2ip and (|22p . So far, the shock description is 
completely general, and depends on the equation of state e(p, v) characterizing the particular 
medium. 

Using Helfand's equation of state (Q and auxiliary relations @ and ([5]), the Hugoniot 
relation (1241) can be used to express in closed form the variation of pressure with specific 
volume across the shock wave, yielding 



7T 



1(1-^) 



a{l-^) 2 -\{l-a) 



(25) 
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FIG. 4. The relation between the shock overpressure and compression ratio represented on the 
shock Hugoniot for a hard disk medium in terms of the initial packing factor rji 

where 7r, a and rji are, respectively, the non-dimensional pressure ratio across the shock 
wave, non-dimensional specific volume ratio across the shock wave, and the initial packing 
fraction of the medium. 



Vi v 2 (Va/m) 
7T = — , a = — , T) X = 

Pi V X Vl 



(26) 



Figure S] displays the shock Hugoniot curves for two initial packing fractions. As can be seen, 
as the packing fraction is increased, a higher pressure is generated for the same compression 
ratio. Alternatively, the same pressure is achieved with less compression in an initially higher 
packed medium. This mimics very well the real properties of shocks in less compressible 
media such as liquids and solids^. The maximum compression that can be achieved across 
a shock wave can be found by letting 7r — > oo in (I25p . This corresponds to letting the 
denominator in (|25|) vanish, yielding a quadratic for which the largest root is the physically 
attainable solution. The maximum compression achievable is thus 



On 



1(1 + 4??i + ^(1 - 8(7?! - l)^; 



6 



(27) 



which is shown graphically in Figure |5j 

It is useful to parametrize the shock jump conditions based on the Mach number of the 
shock wave Mi. From the mass (I21~!) and momentum (|22|) equations and using (111), (151), ([6]) 
and (1261) . we can express the overpressure across the shock wave in terms of the shock Mach 
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FIG. 5. The maximum compression ratio achieved by shock compression in terms of the initial 
packing factor rji 



number and compression ratio across the shock, yielding the so-called Rayleigh line 

7T = 1 + -0(1 - a) 

where ip is simply related to the shock Mach number and the initial value of 7: 

D 2 



(28) 



D 2 

71— = 71 M 

Cl 2 



(29) 



Equating (|25p to (|28|) . we obtain the shock jump conditions in closed form in terms of the 
shock strength parameter if), and the initial packing fraction rji'. 

4 + ip + Aip + Aipr]i 



V2 
Vl 



a 



6ip 



+ 



v/-240(l + tfj)r] 2 + (4 + i) + 4^i) 2 



61/j 



P2 



Z± = 7T= 1 + ^(1 -<T) 
Pi 

D 

T 2 1 - 2i 2 



(30a) 
(30b) 
(30c) 

(30d) 



Figures El [TJ El M show these relations graphically for 771 = 0, r]i — 0.3 and rji = 0.5 in 
terms of the shock Mach number M\. For a given compression ratio a, increasing the initial 
compaction rji gives rise to larger shock Mach numbers. The temperature and pressure 
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FIG. 6. The variation of the compression ratio with the shock Mach number for different initial 
packing factors 
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FIG. 7. The variation of the pressure ratio with the shock Mach number for different initial packing 
factors 

ratios, however, are not very sensitive to changes in the initial compaction levels, as can be 
observed from Figures [7] and 

The parametrization of the shock jump conditions by the shock Mach number, although 
practical, is however not very insightful into the physics of the problem, which correspond to 
the response of a medium to a given prescribed compression. Instead, it is useful to choose 
the square of the piston speed, representing the kinetic energy addition by the piston motion 
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FIG. 8. The variation of the piston to shock speed ratio with the shock Mach number for different 
initial packing factors 




FIG. 9. The variation of the temperature ratio with the shock Mach number for different initial 
packing factors 



as a parameter, given by 



u, 



ei 



(31) 



Expressions ( 1301) can be expressed implicitly in terms of /3 using ( I30c| . which can be rewritten 

as 



P = <!>-. 



:i-a) 



(32) 



(l-^i) 2 

Figures [TUl [TTj [T2l [TBI show the variation of the shock Mach number, specific volume ratio, 
pressure ratio and temperature ratio with j3 varying over more than two orders of magnitude 
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FIG. 10. The variation of the shock Mach number with j3 for different initial packing factors 
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FIG. 11. The variation of the compression ratio with (3 for different initial packing factors 

(shown in logarithmic scale for clarity). For the same piston speed, increasing the initial 
compaction ratio and hence departing from the dilute gas approximation, a weaker shock is 
transmitted with a lower pressure ratio, smaller volumetric compression of the medium and 
lower Mach number. Surprisingly, however, the temperature ratio shown in Fig. [13] for a 
given piston speed is only very weakly affected by the departure from the dilute gas across 
the entire range of weak to strong shocks. Instead, even in the closely packed regimes, the 
temperature ratio is very well approximated by the ideal gas expression, even for the weak 
shocks. To account for this novel finding, the temperature ratio can be written directly from 
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FIG. 12. The variation of the shock overpressure with j3 for different initial packing factors 
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FIG. 13. The variation of the temperature ratio with (3 for different initial packing factors 



the Hugoniot relation (124p which, after some manipulation, can be rewritten as 

T2 _ £2 1 ~ o l u p 2 



(33) 



The second term on the right hand side of fl33l . labeled A, is bounded as u p — > 00, while the 
third term, labeled B, grows without bounds as u p — > 00, hence dominates the temperature 
increase across the shock wave. Formally taking the ratio of these two terms, we obtain 

A _ 2 

B ~ (1 - a) 7 M 2 ( ' 

In the limit of a large shock Mach number, o — > a max (from ( 12?)) ) and hence A/B — > 0. Like- 
wise, this term also vanishes for arbitrary shock Mach numbers when the initial compaction 
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ratio is increased. Indeed, the isentropic exponent increases very rapidly with packing factor, 
as shown in Fig. |2j Under these conditions, the internal energy and temperature, according 
to (1331) . are thus incremented by the kinetic energy of the piston, resulting into an equipar- 
tion of the available energy into the mean motion (mean kinetic energy of the medium) and 
thermal energy. 

To eo 1 u„ 2 , , 

-1 = — £ 1 + — v - 35 
T x ei 2 ei V ; 

This result was also obtained in Ref.— by arbitrarily neglecting the upstream pressure of an 

arbitrary medium in the governing equations. Instead, the condition that the term given in 

(|34|) be small can be considered as the strong shock condition in a hard disk fluid. Generally, 

it was numerically found that the correction for low Mach numbers is very weak, as can be 

seen in Figure [131 suggesting that (|35|) serves as a very good approximation for all initial 

compaction ratios of the medium, ranging from the dilute ideal gas to a solid packed medium 

of hard particles. 



III. COMPARISON WITH MOLECULAR DYNAMIC CALCULATIONS 

In order to validate the results obtained analytically for the shocked state following the 
shock compression by a moving piston, we have performed a series of particle dynamic 
simulations. Owing to the fact that hard particles interact during a vanishingly small time 
compared to their flight time, during which their trajectories are un-affected by any forces, 
their dynamics can be solved in closed forrn^. In between collisions, the particle trajectories 
are known exactly, while the infinitely short collisions only modify the velocity of particles 
upon collision. The collision rules are such that the binary system of colliding particles 
conserves its linear momentum and energy. The sequence of collisions can thus be handled 
quite efficiently in an algorithm evolving the system of many particles collision after collision 
with only very minimal loss in accuracy. This solution strategy is called the Event Driven 
Molecular Dynamics method (EDMD). 

The calculations performed in the present study use the implementation of the EDMD 
method of Poschel and Schwager for two dimensional hard-disk systems^. Consistent with 
the theoretical treatment above, the rotation of the particles is not considered. The code de- 
veloped by Poschel and Schwager was modified to allow for one moving wall, while retaining 
the event-driven nature of the algorithm. The calculations were performed in a rectangular 
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domain of dimensions L x by L y . Particle reflections at all the walls were taken as elastic 
and specular. A total number of 30,000 particles were used. Controlling the dimension of 
the domain set the initial specific volume of the system v\. Adjustment of the radius of the 
particles controlled the initial packing factor of the system rji given by 

NV. Nrtd 2 

The particles were first positioned randomly in the domain, ensuring no overlapping. All the 
particles were assigned equal speeds, with random directions. This initial energy provided 
to the system corresponds to the initial specific energy of the system e\. 

d = \ < V\ > (37) 

which, according to (0Q), sets the initial temperature of the system. The system of particles 
was then let to equilibrate while keeping the number of particles, volume and energy constant 
(i.e., the NVE micro-canonical ensemble), until the distribution of speeds approached the 
Maxwell-Boltzmann distribution. Once the system of particles reached thermal equilibrium, 
the left wall was set in motion with constant velocity u p . An example of the following 
dynamics is shown in Figure O A clearly demarked hydrodynamic shock can be seen 
to propagate across the system of hard disks. In the results reported below, all distances 
have been normalized by the diameter of the hard disks, masses by the atomic mass and 
velocities by the initial speed assigned to the particles U%. Table 1 gives the set of three 
initial conditions investigated, corresponding respectively to a dilute gas (771 = 0.012), dense 
gas (771 = 0.192) and much denser liquid state (r/i = 0.433). Also listed in Table 1 is the 
mean free path, estimated from the kinetic theory of a hard disk medium^ by 

^ __ L x L y 
2V2Nd 

For each packing fraction, the speed of the piston was modified in order to achieve different 
levels of compression. In all simulations, the density and temperature were obtained by 
coarse-grained averages of the state of the gas. At one particular time of the simulation, 
the domain was separated into strips 0.5Ai wide. In each strip, the density (or specific 
volume) was directly obtained by counting particles whose center of mass was within the 
strip. The internal energy, hence temperature, was obtained by taking the average of the 
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FIG. 14. Five consecutive snapshots illustrating the shock wave driven by the moving piston for 
P = 6.25 and rji = 0.192 



kinetic energies of each particle comprised in the strip and subtracting the mean kinetic 
energy in the macroscopic motion, i.e.; 

Uu-(u x }) 2 , N 

T = 2 , / " 39 

k/m 

Because the position of the shock is not known a priori, we have used ensemble averages 
to determine more accurate measurements of the coarse-grained density and temperature. 
This was achieved by repeating the calculation with statistically different initial conditions, 
which were implemented by starting the piston motion at slightly different times after the 
initial equilibration of the medium. An example of the density profile captured during a 
single realization is shown in Figure [151 The ensemble average of 5 of these realizations 
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TABLE I. Parameters for the three fluid regimes investigated via EDMD 
Domain size Packing fraction Mean free path 
L x x L y r]i Ai 

1960 x 1000 0.012 23.1 

490 x 250 0.192 1.44 

327 x 167 0.433 0.642 
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FIG. 15. The density profile captured for a single realization for f3 = 16 and rjx = 0.192 

for density and temperature profiles are shown in Figures [TB] and [T71 As can be seen, the 
density and temperature relax slowly to their equilibrium values behind the shock jump. 
Note that the internal structure of the shock is not within the scope of the present study, 
which deals exclusively with the jump conditions across the shock wave in the thermally 
relaxed media. The results reported below provide the equilibrated values across the shock 
wave. Using the density profiles, we have also determined the position of the shock wave 
at successive times by tracking the position where the density increases by 50% of its total 
increase. This permitted to determine the shock speed. 

Figures dHJ HH] and [2U] show the results obtained in this manner for the specific volume 
decrease across the shock wave, shock wave speed and temperature for different initial com- 
paction ratios 771 and piston speeds u p . Also shown are the analytical predictions using the 
Helfand equation of state detailed in Section 2. In all cases, the agreement was found to 
be excellent, with an error not greater than ~ 3%, even for the large compaction ratio and 

19 



1.5 
I 1 
0.5 


50 100 150 200 

FIG. 16. The average density profile obtained for f3 = 16 and 771 = 0.192 
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FIG. 17. The average temperature profile obtained for /3 = 16 and r/i = 0.192 

strong shocks investigated. Also noteworthy is the fact that the temperature jump across 
the shock wave displayed in Figure [20] shows the predicted invariance with the compaction 
ratio given by (|35|) . 

IV. CONCLUSION 

The compressible dynamics and shock wave propagation in a dense hard disk medium, rel- 
evant to shock propagation in granular hydrodynamics or in liquids, was solved analytically 
using the Helfand equation of state. Physical insight on the role of the medium's initial 
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FIG. 18. Comparison between analytical and calculated compression ratios for different initial 
packing factors and varying shock strength 



-^=0.012 
^=0.192 
- ^=0.433 



FIG. 19. Comparison between analytical and calculated piston speeds for different initial packing 
factors and varying shock strength 

packing on the compressibility of the medium was greatly enhanced by having obtained 
closed form solutions for the isentropic exponent and shock jump conditions. The analytical 
predictions were validated against molecular dynamics calculations using the Event Driven 
Molecular Dynamics technique, where the evolution of the system of colliding elements can 
be obtained analytically. With increasing compaction of the medium, is was shown that the 
medium's compressibility changes substantially, with an isentropic exponent of 7 = 2 in the 
dilute gas phase, and 7 = (9(10) at higher compactions. This was shown to affect signifi- 
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FIG. 20. Comparison between analytical and calculated temperature ratios for different initial 
packing factors and varying shock strength; symbols same as Figure [191 

cantly the shock Hugoniot and shock jump conditions. Parametrization of the shock jump 
relations were obtained using the shock Mach number and piston speed. The important re- 
sult that the temperature in the compressed medium depends to a very good approximation 
only on the square of the piston speed for all compaction levels was shown analytically and 
demonstrated numerically via the molecular dynamic calculations. This important result 
provides a very simple means to estimate the amount of energy injected into a hard particle 
system by surfaces generating shocks. The present study assumed all collisions to be elastic; 
future study will be devoted to studying how the shock hydrodynamics are affected by the 
dissipative nature of granular flows or in reactive flows where collisions may be inelastic. 
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